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Simultaneous electroencephalogram (EEG)-functional Magnetic Resonance Imaging 
(fMRI) recordings have seen growing application in the evaluation of epilepsy, namely in the 
characterization of brain networks related to epileptic activity. In EEG-correlated fMRI stud- 
ies, epileptic events are usually described as boxcar signals based on the timing information 
retrieved from the EEG, and subsequently convolved with a hemodynamic response func- 
tion to model the associated Blood Oxygen Level Dependent (BOLD) changes. Although 
more flexible approaches may allow a higher degree of complexity for the hemodynamics, 
the issue of how to model these dynamics based on the EEG remains an open question. 
In this work, a new methodology for the integration of simultaneous EEG-fMRI data in 
epilepsy is proposed, which incorporates a transfer function from the EEG to the BOLD 
signal. Independent component analysis of the EEG is performed, and a number of met- 
rics expressing different models of the EEG-BOLD transfer function are extracted from the 
resulting time courses. These metrics are then used to predict the fMRI data and to identify 
brain areas associated with the EEG epileptic activity. The methodology was tested on both 
ictal and interictal EEG-fMRI recordings from one patient with a hypothalamic hamartoma. 
When compared to the conventional analysis approach, plausible, consistent, and more sig- 
nificant activations were obtained. Importantly, frequency-weighted EEG metrics yielded 
superior results than those weighted solely on the EEG power, which comes in agreement 
with previous literature. Reproducibility, specificity, and sensitivity should be addressed in 
an extended group of patients in order to further validate the proposed methodology and 
generalize the presented proof of concept. 
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INTRODUCTION 

Over the years, the electroencephalogram (EEG) has been the tool 
of choice for the diagnosis and characterization of epilepsy. With 
the possibility to acquire the EEG simultaneously with functional 
Magnetic Resonance Imaging (fMRI), studies of Blood Oxygen 
Level Dependent (BOLD) signals correlated with epileptic activ- 
ity proliferated (Ives et al., 1993; Hoffmann et al., 2000; Lemieux 
et al., 2001; Salek-Haddadi et al., 2006; LeVan and Gotman, 2009). 
Despite its potential for the localization of epileptogenic brain 
networks in patients with drug-resistant focal epilepsy undergo- 
ing evaluation for surgical treatment, simultaneous EEG-fMRI 
has yet to reach its full potential in clinical practice. One fac- 
tor contributing to this state of affairs is the lack of sensitivity 
in the identification of hemodynamic changes associated with 
the EEG epileptiform discharges in a significant number of stud- 
ies (Aghakhani et al., 2006; Salek-Haddadi et al., 2006; Gotman, 
2008; Grouiller et al., 201 1). Although technical difficulties related 
with data acquisition and artifact correction may in part explain 
such results, probably most important are the limitations related 
with the remaining conceptual and methodological challenges 
associated with the integration of the two types of signals. 



Although a great amount of both experimental and theoretical 
work has been dedicated to the clarification of the relationship 
between neural activity and associated hemodynamic changes, 
neurovascular coupling mechanisms remain an active area of 
research (Rosa et al., 2010a). The most consensual evidence comes 
from the recording of electrical activity using micro-electrodes 
implanted in the cortex of (non-human) animals, simultane- 
ously with fMRI, indicating that the BOLD signal reflects mostly 
slow, post-synaptic input activity measured by local field poten- 
tials (LFPs), rather than fast, spiking output activity measured 
by single/multi-unit activity (S/MUA; Logothetis et al., 2001). In 
humans, a growing number of simultaneous EEG-fMRI studies 
on healthy subjects as well as epilepsy patients have now been 
reported (Goldman et al., 2002; Laufs et al., 2003, 2006; Moos- 
mann et al., 2003; de Munck et al., 2009), and biophysical models 
of the neurovascular coupling have been proposed (Riera et al., 
2006, 2007). Overall, reports in the literature do not provide a 
clear picture of the link between EEG and BOLD signals. In par- 
ticular, contradictory results have been presented regarding the 
dependency of BOLD changes on the EEG power and spectral 
profiles. These include, for example, positive and negative BOLD 
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correlations with specific frequency band power changes in the 
human EEG (de Munck et al., 2009), BOLD decoupling from LFP 
power in mice (Ekstrom, 2010), and negative BOLD associated 
with large increases in LFP and MUA during seizures also in mice 
(Schridde et al, 2008). Rosa et al. (2010b) addressed this topic 
by comparing different models of the transfer function between 
EEG and BOLD signals, in the prediction of fMRI data, in a visual 
stimulation experiment with human healthy subjects. The mod- 
els explored included the EEG total power (TP; Wan et al, 2006); 
linear combinations of the power from different frequency bands 
(Goense and Logothetis, 2008); and several variations of a heuris- 
tic model proposed by Kilner et al. (2005) in which BOLD changes 
are assumed to be proportional to the root mean square frequency 
(RMSF) of the EEG spectrum. The results obtained showed a clear 
superiority of the RMSF metrics in predicting the BOLD signal, 
when compared to power-weighted metrics. 

In most epilepsy EEG-fMRI studies, the goal is to find brain net- 
works exhibiting hemodynamic changes associated with interictal 
and/or ictal activity, which are expected to be correlated with the 
epileptogenic areas (Hoffmann et al., 2000; Salek-Haddadi et al., 
2006; Marques et al., 2009). Both ictal and interictal events are 
identified on the EEG trace and are then used to define regressors 
of interest in a general linear model (GLM) analysis of the fMRI 
data. Interictal spikes are usually described as stick functions and 
ictal activity as boxcar signals between seizure onset and offset, 
eventually sub-divided into up to three phases: early ictal EEG, 
clinical onset, and late ictal EEG. Both types of events are then 
convolved with a model of the hemodynamic response function 
(HRF; Tyvaert et al, 2008; Moeller et al, 2010; Thornton et al, 
2010). Independent component analysis (ICA) of the fMRI data 
has also been performed in order to identify interictal/ictal BOLD 
patterns without resorting to the EEG (LeVan et al., 2010; Thorn- 
ton et al., 2010). When the EEG accurately reflected the seizure 
onset, the GLM approach yielded activations concordant with the 
ictal onset zone; otherwise ICA still gave valuable insight on the 
ictal hemodynamics. Importantly, the question of whether or not 
the neurovascular coupling is preserved from healthy to disease 
conditions has also been investigated, by allowing variations in the 
HRF (Grouiller et al, 2010). It is in principle possible to achieve 
a higher degree of complexity for the epileptic activity hemody- 
namics by using more flexible HRFs, or completely model-free 
approaches such as ICA. However, such approaches do not address 
the issue of how to model these dynamics based on the information 
available in the EEG data (Lemieux, 2008). 

In this paper, we address the issue of modeling the BOLD 
dynamics associated with epileptic EEG activity by extracting met- 
rics from the EEG spectrum expressing different models of the 
transfer function between neuronal and hemodynamic signals. A 
methodology is proposed for the analysis of EEG-fMRI data in 
epilepsy, consisting of ICA decomposition of the EEG followed 
by component selection based the reproducibility across different 
acquisition runs, Morlet wavelet spectral analysis, and EEG metric 
extraction. The resulting time courses are convolved with a canon- 
ical HRF and used as regressors of interest in a GLM analysis of 
the fMRI data. The proposed methodology is applied to the study 
of a patient with epilepsy associated with a hypothalamic hamar- 
toma. The different EEG-fMRI transfer functions are compared 



with each other, as well as with a conventional GLM methodology 
based on the identification of ictal and interictal events on the EEG 
by the neurophysiologist, and also with an fMRI-ICA approach. 

MATERIALS AND METHODS 
PATIENT CHARACTERIZATION 

We focused on the simultaneous EEG-fMRI data recorded from a 
2-year-old patient with a giant hypothalamic hamartoma suffer- 
ing from gelastic epilepsy, as part of the pre-surgical evaluation 
under the Program for Epilepsy Surgery of the Hospital Center of 
West Lisbon. This case study has been previously described (Leal 
et al., 2009). This patient was studied in a run of 30 patients, 
from which only five had ictal events during the scanning ses- 
sions. From these five, only this one patient had an EEG trace clear 
from movement related MR artifacts triggered by the beginning 
of the seizures and therefore presented sufficient data quality for 
subsequent EEG quantification and was selected for this study. 

Seizures occurred more than 50 times per day and typically 
lasted for 20-30 s, involving almost exclusively the left hemisphere. 
The clinical manifestations of the seizures consisted of slowing 
of motor activity, variable interruption of consciousness, eye- 
lid rhythmic movements with bilateral nystagmus to the right, 
and occasionally gelastic laughter. After the acquisition of the 
EEG-fMRI data, the patient underwent a two-stage hamartoma 
disconnection surgery, 1 year after which the seizures were reduced 
to 1-3 per day. 

The EEG interictal activity demonstrated a persistent slow- 
wave abnormality over the left temporal-occipital area, associated 
with abundant spike activity with occasional contralateral propa- 
gation. Abundant spike activity also occurred over the left hemi- 
sphere frontal lobe. Topological analysis of the interictal spikes 
presented a spatial stationarity for the posterior spikes, whereas the 
frontal ones changed significantly in configuration from an occip- 
ital dipolar potential at spike onset to a dipolar frontal potential at 
spike peak. The ictal EEG pattern was very monotonous and con- 
sisted of early diffuse desynchronization, followed by the build-up 
of spike activity over the left occipital and temporal areas and, in 
the later stages of the seizure, over the frontal area. Occasionally 
secondary propagation of spike activity to the right temporal areas 
occurred. 

SIMULTANEOUS EEG-fMRI ACQUISITION 

The EEG was recorded using an MR-compatible 37-channel sys- 
tem (Maglink, Neuroscan, Charlotte, NC, USA) with two ECG 
channels, using a sampling rate of 1000 Hz with a bandwidth DC- 
250 Hz and reference at electrode FCz. The imaging was performed 
on a 1.5 T MRI scanner (GE Cvi/NVi). Six fMRI runs were col- 
lected using a gradient-echo echoplanar imaging (EPI) sequence, 
with TR = 2.275 s, 3.75 mm x 3.75 mm x 5.00 mm voxel size and 
a total of 154 volumes (the first four were subsequently rejected in 
order to discard Ti relaxation unstationarities). A Ti -weighted 
structural image was also acquired using a 3D spoiled gradi- 
ent recovery (SPGR) sequence, with 0.94 mm x 0.94 mm in-plane 
resolution and 0.6 mm slice thickness. During the scanning ses- 
sion, the patient was administered light anesthesia with Sevoflu- 
rane at 1% (Abbott Laboratories, Abbot Park, IL, USA), through 
mask, as established by the MRI protocol for small children and 
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uncooperative patients. A 5min EEG recording was also per- 
formed outside the scanner and before anesthesia, in order to allow 
for cross-validation of the simultaneous EEG-fMRI recordings in 
terms of the presence of MR related artifacts and the effects of 
anesthesia. Periods of interictal and ictal activity were identified 
on the artifact-corrected EEG by the neurophysiologist. For runs 
2 and 3, two and five ictal events occurred, respectively. For the 
remainder runs, only interictal spikes were detected on the EEG. 
The visual inspection of the EEG inside and outside the scanner 
revealed a slight increase in both beta and theta background activ- 
ity under anesthesia, but there was no significant change in the 
morphology of interictal spikes. 

EEG ANALYSIS 
Pre-processing 

The EEG pre-processing was executed using the EEGLAB toolbox 
(Delorme and Makeig, 2004). Firstly, the EEG traces were visu- 
ally inspected for the presence and consequent rejection of bad 
channels, associated with poor contacts. A 2 Hz high-pass filter 
was then applied so as to remove baseline drifts from the signal. 
The times of occurrence of the gradient artifacts associated with 
the acquisition of each fMRI slice were automatically identified on 
the EEG signals, by software developed in-house. The fMRI gra- 
dient artifact correction algorithm, FMRIB's FASTR (Niazy et al., 
2005), was then applied using the default parameters. For the pulse 
artifact removal, FMRIB's QRS complex identification algorithm 
(Niazy et al, 2005) was first applied to the ECG channels. An opti- 
mal basis set of three principal components was then employed for 
pulse artifact removal of the data, after low-pass filtering at 45 Hz 
and down-sampling to 100 Hz to improve manageability. 

ICA decomposition 

In an attempt to separate out the activity of interest in the EEG, 
the pre-processed data were decomposed by ICA using the info- 
max algorithm as implemented in EEGLAB (Delorme and Makeig, 
2004). The reference channel was arbitrarily kept as the one cho- 
sen by the electrophysiologist during the acquisition. Although the 
referencing method for the EEG channels does not affect the final 
IC time courses, because the reference channel is linearly separable 
from the data, it will affect the IC's scalp topographies. Neverthe- 
less, the reference channel was kept the same (FCz) throughout 
the acquisition of all six runs, so comparisons between the scalp 
topographies of components obtained in different runs are still 
possible. 

A reproducibility analysis of the ICs was performed in order 
to identify the associated topographies that were consistent across 
the six acquisition runs. IC groups were built by fixing each IC 
of each run and selecting, for each of the other runs, the IC that 
was the most spatially correlated with it. A total of 6 runs x 25 
ICs= 150 IC groups were hence generated, each composed of a 
string of six ICs. An IC group was considered to be consistent if 
the same associated string was generated by one IC of each run, and 
therefore was repeated six times in the total set of 150 IC groups. 
The topographies and time courses of the consistent IC groups 
were then visually inspected for the identification and consequent 
rejection of artifact related ICs, dominated by residual gradient 
artifacts, bad channels, or eye blink/movement artifacts. 



EEG metrics 

The spectral profiles of the selected ICs were obtained by time- 
frequency analysis through convolution of the respective time 
courses with Morlet wavelets. The subset of metrics found to 
be more relevant in Rosa et al. (2011) were applied here: mean 
frequency (MF), RMSF, un-normalized mean frequency (uMF) 
un-normalized mean square frequency (uRMSF), and TR 

The difference between frequency averaging measures (RMSF 
or MF) was found to have a negligible effect on BOLD signal pre- 
diction; hence the results emerging from the MF and uMF metrics 
will be omitted, as they were not significantly different from those 
of the RMSF and uRMSF metrics, respectively. 

fMRI ANALYSIS 

The fMRI data were analyzed using FSL 1 , including: (1) pre- 
processing; (2) GLM; and (3) ICA. 

Pre-processing consisted on: motion correction using 
MCFLIRT (Jenkinson et al., 2002); slice timing correction using 
(Hanning-windowed) sine interpolation to shift each time-series 
by an appropriate fraction of a TR relative to the middle of the TR 
period; brain extraction using BET (Smith, 2002); temporal high- 
pass filtering rejecting periods above 100 s; and spatial Gaussian 
filtering with FWHM = 8 mm. 

Two GLM analyses were specified in order to identify BOLD sig- 
nal changes associated with: ( 1 ) the EEG metrics for each consistent 
IC group (as defined in the previous section); and (2) the epileptic 
activity identified by the neurophysiologist on the pre-processed 
EEG traces, where boxcar signals were used to describe the periods 
of ictal activity and stick functions were used to describe interictal 
spikes (from now on referred to as EA). Each of these variables 
of interest (EEG metrics and EA) was convolved with a canonical 
HRF (Friston et al, 1998), and the time and dispersion deriva- 
tives were also included to account for some degree of variability 
in the HRF shape across the patient's brain. This resulted in a 
set of three regressors for each variable. The final GLM regres- 
sors were obtained by re-sampling the resulting time courses to 
match the middle of the acquisition time period of each fMRI 
volume. Six motion parameters were also included as regressors 
of no interest, in order to account for residual motion-related sig- 
nal jitter not removed by the motion correction procedure. The 
GLM's were fitted to the data using the FILM algorithm (Woolrich 
et al, 2001) F tests were then applied to each estimated para- 
meter contrast, resulting in Z (Gaussianized F) statistic maps. 
These were thresholded using a clustering procedure, whereby 
each cluster is determined by a voxel Z > 2.3 and a (corrected) 
cluster significance threshold p = 0.05. 

An ICA of each fMRI run was also performed using Probabilis- 
tic ICA as implemented in MELODIC (Multivariate Exploratory 
Linear Decomposition into Independent Components) Version 
3.09, part of FSL (FMRIB's Software Library (see text footnote 
1); Beckmann and Smith, 2004). Pre-processing included voxel- 
wise de-meaning of the data and normalization of the voxel-wise 
variance. Estimated component maps were divided by the stan- 
dard deviation of the residual noise and thresholded by fitting a 
mixture model to the histogram of intensity values. 



1 www.fmrib.ox.ac.uk/fsl 
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MODEL COMPARISONS 

The comparison between EEG metrics was performed through 
GLM analysis followed by inferences based on F tests, for each 
consistent IC group of interest and for each acquisition run. We 
first considered a single GLM comprising the three EEG met- 
rics (TP, RMSF, and uRMSF) with three regressors of interest 
each (canonical HRF, time, and dispersion derivatives), totaling 
nine regressors of interest, in one three-way comparison. Three 
GLM's contrasting pairs of EEG metrics (TP vs. RMSF, TP vs. 
uRMSF, and RMSF vs. uRMSF) were also considered, totaling 
six regressors of interest each, in three two-way comparisons. 
F tests were computed for each set of three regressors regard- 
ing each metric, as well as for the whole set of metrics. The 
resulting Z statistical maps were inspected for their volume, i.e., 
the total number of voxels exhibiting significant EEG metric- 
related BOLD changes. This was used as a quantitative measure 
of the performance of each EEG metric in terms of BOLD pre- 
diction, and the analysis was repeated for each consistent IC 
group. 

The comparison between consistent IC groups of interest was 
performed in a similar way, for each acquisition run, applying the 
selected EEG metric. Here, we considered a single GLM compris- 
ing the three consistent IC groups of interest (I, II, and V, as will 
be shown in the Results) with three regressors of interest each 
(canonical HRF, time, and dispersion derivatives), totaling nine 
regressors of interest, in one three-way comparison. 



In order to assess the plausibility and consistency of the EEG 
metric-derived BOLD maps, the EEG metric/IC group combina- 
tion yielding the largest maps (and hence best at predicting BOLD 
signal changes) was selected for subsequent comparison with the 
GLM analysis using the EEG EA regressors and also with the ICA 
analyses. The consistency between each two regressors was mea- 
sured as their temporal correlation. The consistency between two 
maps was measured as their spatial overlap, i.e., the ratio of the 
volume of the map intersection with the volume of the map union. 

RESULTS 

In this section, the results obtained through the proposed method- 
ology will be presented. 

EEG ANALYSIS 

An exerpt of the EEG trace obtained after pre-processing is pre- 
sented in Figure 1. The data appears to be clear from MR related 
artifacts and the ictal activity can be clearly identified. 

The consistent IC groups, obtained as a function of the repro- 
ducibility of the associated IC topographies across runs, are pre- 
sented in Figure 2. The IC groups I, II, and V were considered 
artifact free and were kept for further analysis, while the remainder 
were rejected. The IC groups I and II exhibit clear dipolar config- 
urations in the left hemisphere, compatible with the frontal and 
occipital/parietal patterns of the patient's interictal and ictal EEG. 
Given its predominantly frontal topography, particular attention 
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FIGURE 1 | Example of the EEG trace obtained after pre-processing the data for run 3. An ictal event occurs within the time window starting at 182 s and 
ending at 193 s, as indicated. The data appear to be clear of MR related artifacts. 
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FIGURE 2 | Consistent IC groups I— VII: the associated IC topographies for each run are shown (color scale indicates relative weight of each EEG 
channel), together with the respective IC number. The three IC groups highlighted in blue boxes were considered not artifactual and hence of interest (left is 
left, right is right, up is anterior, down is posterior). 



was given to IC group II in order to confirm that ocular move- 
ment artifacts did not dominate the IC time course. Topogra- 
phy group V exhibits a more diffuse configuration difficult to 
interpret. 

The spectral profiles obtained with Morlet wavelet decomposi- 
tion for one IC (IC10) in run 3, as well as the EEG channel with 
highest absolute weight for this component (P3), are shown in 
Figure 3. Spectral changes associated with the ictal events iden- 
tified by the neurophysiologist are visible in both spectrograms; 
however, these changes are clearer in the spectrogram obtained 
for the IC in comparison to the spectrogram of channel P3, or 
any other single EEG channel (data not shown). This observation 
suggests that ICA decomposition is capable of separating out 



the ictal activity spread over several EEG channels into a limited 
number of components. 

fMRI ANALYSIS 

The BOLD statistical maps obtained using each EEG metric and 
consistent IC group were generally consistent with each other, 
across runs and also with the patient's seizure semiology, but 
differed considerably in terms of the number of voxels showing 
statistically significant EEG metric-related BOLD changes. Firstly, 
the results of the comparison between the fMRI analysis using the 
different EEG metrics and consistent IC groups will be presented. 
A single metric and IC group will then be selected for comparison 
with the EA GLM and fMRI-ICA analyses. 
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FIGURE 3 | Spectrograms of the time courses of channel P3 (top) and IC 10 (group II; middle) with corresponding EEG metrics (bottom), for the EEG of 

run 3. The P3 channel was the one that contributed the most to IC 10. The boxes in red/black represent the periods identified as ictal events. 



EEG metric comparisons 

The three-way analysis did not show meaningful activations to 
allow for the comparison of the individual EEG metrics; however, 
the two-way comparisons yielded a clear superiority of the RMSF 
metric when compared to the TP and uRMSF metrics, for every 
IC group and run, in terms of the number of voxels showing sta- 
tistically significant EEG metric-related BOLD changes, as shown 
in Figure 4. 

The three-way comparison regarding the IC group showed a 
significant superiority of IC group II, in terms of the number 



of voxels exhibiting statistically significant EEG metric-related 
BOLD changes, as shown in Figure 5. For five out of the six runs, 
the statistical maps yielded by IC group II had larger volumes 
of significant voxels than those yielded by other ICs. The EEG 
RMSF metric of IC group II was found to be the best at predict- 
ing BOLD changes and it will now be compared with the results 
of the EA GLM analysis as well as with those of the fMRI-ICA 
approach. 

The BOLD statistical maps obtained using the RMSF metric 
for IC groups I and V are shown in Figures 6 and 7. 
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FIGURE 4 | Number of voxels exhibiting statistically significant EEG 
metric-related BOLD changes, for the two-way comparisons RMSF vs. 
TP (top) and RMSF vs. uRMSF (bottom), for each IC group (II, V, and I) 
and each run (1-6). 



EEG metric vs. EA GLM analysis 

The comparison between the fMRI results obtained using the EEG 
RMSF metric and the EA GLM analysis, for IC group II, are sum- 
marized in Figure 8. For each run, the regressors associated with 
the canonical HRF are plotted alongside with the Z (Gaussianized 
F) statistic maps obtained with both methods. The temporal cor- 
relation between the regressors and the spatial overlap between the 
respective maps are also presented. The consistency between the 
two methodologies is generally low for all runs (correlation < 0.3 
and overlap < 3), with only run 3 exhibiting an overlap between 
maps above 10%. Interestingly, this is also the run during which 
more ictal events occurred. 

In general, the maps obtained using the EEG metric approach 
were consistent across runs and also with the patient's seizure semi- 
ology. In fact, clusters were found on the left occipital/parietal and 
left frontal lobes, consistently with the observation on the EEG 
of spike buildups in left occipital/parietal channels followed by 
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FIGURE 5 | Number of voxels exhibiting statistically significant 
RMSF-related BOLD changes, for the three-way analysis of IC groups I, 
II, and V, for each run (1-6). 



the left frontal channels, after a generalized desynchronization. 
In run 3, the EA approach was also capable of identifying these 
brain areas, but the corresponding statistical maps were more sig- 
nificant and extensive for the EEG metric approach. Furthermore, 
with the EEG metric approach, clusters involving left thalamic, left 
hippocampal, and left frontal ventral areas, as well as the hamar- 
toma itself, were also found, which are generally consistent with 
the brain network associated with seizures originating in a hypo- 
thalamic hamartoma. For the remaining runs, in contrast with the 
EEG metric approach, the fMRI statistical maps obtained with the 
EA approach were in general not consistent across runs nor with 
the expected epileptic network. 

EE G metric vs. fMRI-ICA approach 

The results regarding the EEG metric GLM analysis and the fMRI- 
ICA decomposition of all the runs, for IC group II, are presented 
in Figure 9. The fMRI-ICs presented are those that yielded the 
highest spatial map overlap with the results of the correspond- 
ing EEG metric analysis. Map spatial overlaps and time course 
temporal correlations between EEG metric GLM and fMRI-ICA 
approaches were generally higher than those found between EEG 
metric GLM and EA GLM approaches (Figure 8). Moreover, in 
contrast with the EA approach, with fMRI-ICA, consistency across 
runs was also observed for the brain areas expected in the patient's 
epileptic network. Interestingly, runs where only interictal activity 
was recorded (1, 4, 5, and 6) yielded maps consistent with runs 
with ictal activity (2 and 3). 

DISCUSSION 

We proposed a new methodology for BOLD signal prediction in 
EEG-correlated fMRI studies in epilepsy, by incorporating a model 
of the EEG-BOLD transfer function. Specifically, independent 
components of the EEG associated with consistent topographies 
were translated into BOLD signal predictions by a set of model- 
based metrics. Interestingly, we found that increases in the MF of 
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FIGURE 6 | Time courses (left) and Z statistic maps (right) of the RMSF metric for IC group I. 



the EEG were better than its power at predicting BOLD increases, 
in support of the heuristic proposed in (Kilner et al, 2005) and 
in agreement with the results obtained in a visual stimulation 
experiment with healthy subjects (Rosa et al., 2010b). Moreover, 
such EEG-based metrics were found to improve detection sensi- 
tivity compared with conventional approaches to EEG-fMRI data 
analysis in epilepsy. 

The heuristic proposed by Kilner and colleagues puts forward 
the MF (in an RMS sense) of the EEG spectrum as a surrogate of 
the neuronal activity eliciting the BOLD signal, following a broad 
physiological inspiration: the BOLD eliciting signal is assumed to 
be proportional to the electric work dissipated by the ionic cur- 
rents across the cell membranes; this can in turn be shown to be 
proportional to the RMSF of the LFP/EEG spectrum if the covari- 
ance of the membrane potentials of the cells are assumed constant 
(Kilner et al, 2005). Although more detailed models of EEG-fMRI 



integration have been presented in the literature (Riera et al, 2006, 
2007), the heuristic benefits from its simplicity in application and 
interpretability. The dependency of the BOLD signal on the EEG 
spectral profile has often been experimentally reported in stud- 
ies of spontaneous or evoked fluctuations of brain rhythms, and 
the results are most frequently found to be in concordance with 
Kilner 's heuristic predictions (Goldman et al., 2002; Goncalves 
et al, 2006; Laufs et al, 2006; de Munck et al, 2009; Michels 
et al, 2010; Zumer et al, 2010; Khursheed et al, 2011). Rosa 
et al. (2010b) have explicitly employed the heuristic to analyze 
EEG-fMRI data collected during a visual flicker stimulation task 
in healthy subjects, and showed that BOLD signal decreases were 
indeed associated with changes in the EEG spectral profile, namely 
its RMSF, which did not arise from power changes in one spe- 
cific frequency band. In epilepsy, low- frequency slow- wave activity 
increases have been shown to be associated with BOLD decreases 
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FIGURE 7 | Time courses (left) and Z statistic maps (right) of the RMSF metric for IC group V. 



(Archer et al., 2003) while high-frequency spike and wave dis- 
charges have been shown to be associated with BOLD increases 
(Krakow et al, 2001; Hamandi et al, 2004). Although these find- 
ing are in agreement with the heuristic, our study is the first one 
to test it explicitly on EEG-fMRI data of epileptic activity. 

The EEC metric-related BOLD change maps were consistent 
with the ones obtained using the epileptic activity regressors 
defined by the neurophysiologist, whenever ictal activity was iden- 
tified on the EEC For the runs on which no ictal events were 
detected on the EEC, the proposed methodology was able to iden- 
tify the same brain network that was involved in the ictal runs. This 
was achieved by regressors based on the same IC scalp topography 
as those from the ictal runs, suggesting the presence of under- 
lying epileptic activity in the identified epileptic network, which 
only occasionally manifests itself with an ictal character. In con- 
trast with the proposed EEC metric based approach, however, the 



conventional analysis of the interictal runs yielded inconsistent 
or no results for the same statistical significance threshold. These 
findings suggest that the interictal events detected on the EEC may 
not fully reflect the activity of the underlying brain network, while 
the selected EEC metrics may be more powerful in depicting it. 
The same network was also identified by fMRI-ICA on both ictal 
and interictal runs, which further supports this idea. However, the 
fully data-driven method of fMRI-ICA lacks an explanatory model 
for the data, in contrast with the proposed methodology, which is 
based on modeling the EEG-BOLD transfer function. 

For the purpose of verifying the specific epileptic character of 
the identified brain networks, their consistency across runs and 
their plausibility as the epileptic brain network underlying the 
patient's seizure semiology were considered. The overlap of the 
corresponding BOLD change maps across runs was very high, 
both for the EEG metrics GLM and the fMRI-ICA approaches. 
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Moreover, these maps were in good agreement with the brain net- 
work known to be involved in the epileptic activity of this patient, 
which includes the hamartoma, as well as left hemisphere hypo- 
thalamus, hippocampus, parietal-occipital lobe, cingulate gyrus, 
and dorsal-lateral frontal lobe (Leal et al., 2009). Clusters in the 
left parietal-occipital and frontal lobes are consistent with the 
observation on the patient's EEG of spike buildups in left parietal- 
occipital channels followed by the left frontal channels, after a 
generalized desynchronization. Clusters involving left thalamus 
and hippocampus, as well as the hamartoma itself, are generally 
consistent with the brain network associated with seizures orig- 
inating in a hypothalamic hamartoma (Leal et al., 2003). Future 
work should be carried out in order to further validate the net- 
works identified by our methodology. On a first approach, the 



topographies of the selected ICs could be compared with those 
obtained by ICA decomposition of the EEG performed outside the 
MRI scanner. Ultimately, intra-cranial EEG recordings (unavail- 
able for this patient) must be used in order to achieve a conclusive 
validation. 

An ICA decomposition of the EEG was used here with the pur- 
pose of separating out the activity of interest into a univariate 
timecourse, which was expected to exhibit a consistent and mean- 
ingful topography. ICA is a popular technique for the removal of 
muscular activity or eye movement artifacts in EEG data process- 
ing (Vigario, 1997). Because of the large amplitude of interictal 
epileptic activity and the fact that its sources can generally be 
assumed to be spatially stationary, ICA has also proved to be use- 
ful in the separation and identification of such activity (Kobayashi 
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et al, 1999; Urrestarazu et al., 2006; Marques et al, 2009; Formag- 
gio et al., 201 1). However, ICA decomposition of EEG ictal events 
that involve spatial propagation may be questionable in the sense 
that the spatial stationarity assumption of the EEG sources is not 
verified. This is in fact the case in our study, where ICA was applied 
to the EEG recorded during seizures exhibiting a spatial propaga- 
tion pattern associated (Leal et al., 2009). The results obtained may 
reflect this issue to some extent, as no single IC isolated per se all 
of the seizure dynamics. Nevertheless, ICA was useful for the sep- 
aration of local approximately stationary activity, giving the ICs 
higher signal-to-noise ratio for the signals of interest, by separat- 
ing them from activity in neighboring brain regions and also from 
residual artifacts not fully corrected by the EEG pre-processing as 
observable in the spectral analyses and scalp topographies. 



Other approaches for the extraction of a univariate EEG time 
course, representative of the epileptic activity, have been proposed 
in the literature, namely continuous Electrical Source Imaging 
(Vulliemoz et al., 2010) and weighted averaging of selected ICs 
(Formaggio et al., 201 1). The former corresponds to the projection 
of the recorded EEG data into the space of an EEG source estima- 
tion, which provides a more informed way of obtaining the EEG 
source signal than ICA. However, this technique is more prone 
to include artifacts in the resultant time courses when compared 
to ICA, as the latter automatically rejects channels contaminated 
with relevant artifacts. Regarding the averaging of selected EEG 
IC time courses, there is no straightforward way to compute 
averaging weights. Furthermore, when averaging ICs with dis- 
crepant topographies, one incurs in the risk of mixing sources 
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with significantly different dynamics and loosing the meaning of 
the ICA source separation. 

A limitation of the proposed approach is the bias of the EEG 
measures toward superficial cortical activity, which possibly pre- 
cludes the identification of hemodynamic changes associated with 
deep brain activity. This limitation is however common to all EEG- 
fMRI "integration through prediction" approaches. Nonetheless, 
the particular syndrome of epilepsy associated with hypothala- 
mic hamartomas is rather well described in the literature, with 
clear evidence for the hypothalamic hamartoma being the seizure 
focus, and this region was in fact found in our work in five out of 
six EEG-fMRI datasets. Although the proposed EEG-BOLD trans- 
fer functions are not specific to epileptic activity, this specificity 
is achieved in the presented methodology by selecting the EEG 
topography used to extract the metric based on an ICA proce- 
dure. Nevertheless, a possible limitation of the proposed transfer 
functions is that they are not specific to the start of the seizure, 
and hence to the epileptogenic focus. However, our aim was the 
description of the full epileptic network, leaving the specific iden- 
tification of the seizure focus and seizure propagation dynamics 
for other, related lines of research (Murta et al., 2012). 

The study presented here focused on data from a single patient 
with the aim of providing a proof of concept for the poten- 
tial usefulness of the proposed methodology for the identifica- 
tion of epileptic networks. The choice of this patient was based 
on the quality of the EEG data that could be achieved due to 
the absence of movement associated with the beginning of the 



seizures. Moreover, epilepsy cases associated with hypothalamic 
hamartomas are relatively stereotypical in terms of the electro- 
physiological patterns of seizure propagation, which makes the 
interpretation of the results relatively more robust in comparison 
to other types of epilepsy. In general, the clinical utility of EEG- 
fMRI is yet to be established. Nevertheless, in this case study, the 
results of the EEG-fMRI investigation indicate possible alterna- 
tive treatment approaches involving the surgical interruption of 
the seizure propagation pathways (Leal et al., 2009; Murta et al, 
2012). The proposed methodology should now be applied to an 
extended group of patients, in order to generalize the proof of 
concept presented here and to further validate it. Reproducibility, 
specificity, and sensitivity should then also be addressed. 

In conclusion, we presented a new approach for EEG-fMRI 
integration in the field of epilepsy, which incorporates and tests 
different models of the transfer function between EEG and BOLD 
signals, hence allowing better predictions of the hemodynamic 
changes associated with epileptic activity. This work therefore pro- 
vides a contribution to our understanding of the link between EEG 
and BOLD signals as well as for improving the yield of EEG-fMRI 
studies in epilepsy. 
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